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ABSTRACT 

We present a new approach to studying the evolution of massive black hole binaries in a stellar 
environment. By imposing conservation of total energy and angular momentum in scattering 
experiments, we find the dissipation forces that are exerted on the black holes by the stars, 
and thus obtain the decaying path of the binary from the classical dynamical friction regime 
down to subparsec scales. Our scheme lies between scattering experiments and A^-body sim- 
ulations. While still resolving collisions between stars and black holes, it is fast enough and 
allows to use a large enough number of particles to reach a smooth and convergent result. We 
studied both an equal mass and a 10:1 mass ratio binaries under various initial conditions. We 
show that while an equal mass binary stalls at a nearly circular orbit, a runaway growth of 
eccentricity occurs in the unequal mass case. This effect reduces the timescale for black hole 
coalescence through gravitational radiation to well below the Hubble time, even in spherical 
and gasless systems formed by dry mergers. 
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1 INTRODUCTION 

Supermassive black hole (BH) binaries are formed as a result of 
galactic mergers. The two BHs sink towards the bottom of the po- 
tential well and form a bound pair at the centre of the stellar distri- 
bution. According to the classical picture (Begelman, Blandford & 
Rees 1980) the binary hardens until the loss cone is depleted (i.e. 
depleted of stars on low angular momentum orbits that pass close 
enough to the binary for a significant interaction), and stalls at a 
separation too large for the emission of gravitational radiation to 
cause rapid inspiral and coalescence (the 'final parsec problem'). 
This picture has been tested extensively using N-body simulations 
and scattering experiments. Early numerical studies concentrated 
on the simple case of an equal mass binary on a circular orbit 
within a spherically symmetric stellar distribution. Those studies 
(e.g. Makino 1997; Quinlan & Hernquist 1997; Milosavljevic & 
Merritt 2001; 2003), with a relatively small number of particles, 
were somewhat inconsistent with each other and could not confirm 
loss cone depletion due to spurious numerical relaxation over long 
time scales. 

In recent years, more detailed investigations have been per- 
formed. Apart from the increased number of particles and longer 
integration times, various complications have been considered such 
as triaxiality and rotation (e.g. Berczik et al. 2006; Khan, Just & 
Merritt 2011), post Newtonian corrections (e.g. Berentzen et al. 
2009) massive perturbers (e.g. Perets & Alexander 2008) and gas 
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discs (e.g. Cuadra et al. 2009). These studies found that binary co- 
alescence in less than the Hubble time is feasible even if only one 
of these factors is present. Thus, the final parsec is potentially not a 
problem in a realistic merger remnant. 

More recent work (e.g. Iwasawa et al. 2011; Sesana, Gualan- 
dris & Dotti 2011) has focused on the evolution of the binary 
eccentricity rather than just the binary separation. For high mass 
ratios (~ 100:1), these authors showed that even in a spherically 
symmetric environment, eccentricity can increase to almost unity. 
Very high eccentricity means that in pericentre the two BHs can 
be close enough together that gravitational wave emission becomes 
efficient. This enhanced energy loss to gravitational radiation (as 
compared to the circular case with the same semi-major axis) can 
reduce coalescence timescale to well below the Hubble time. 

Direct W-body simulations are the most accurate way to study 
binary BH (BBH) evolution, but they are computationally expen- 
sive. It is therefore difficult to perform diverse enough tests to cover 
the problem's parameter space. By compromising for an unrealisti- 
cally small number of particles, one introduces spurious relaxation. 
In spherical galaxies, this process drives unrealistic loss cone re- 
population (Yu 2002; Milosavljevic & Merritt 2003) and causes the 
hardening rate to be highly N-dependent (Makino & Funato 2004; 
Berczik, Merritt, & Spurzem 2005). Berczik et al. (2006) followed 
the evolution of a BBH in triaxial and rotating galaxy models and 
found that the hardening rate was iV-independent, implying a colli- 
sionless mode of loss cone repopulation. 

Our previous work dealt with a possible stellar kinematical sig- 
nature of stars around a BBH. In Meiron & Laor (2010) we pro- 
duced kinematical maps (projected maps of mean stellar veloc- 
ity, velocity dispersion and higher velocity moments); this was 
achieved by scattering stars on a BBH on a fixed circular orbit. 
To produce more realistic kinematical maps, we looked for a way 
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to calculate a BBH path fast and accurately so it could be used for 
following scattering experiments. Good kinematical maps require a 
very large number (N < 10 8 ) of stars on small scales, otherwise the 
high moments of the line of sight velocity distribution are poorly 
resolved. A'-body simulations cannot be made yet with such a large 
number of particles on small enough scales, but this is not a prob- 
lem for scattering experiments which are performed on a precal- 
culated BBH path. Since we were already working with scattering 
experiments, we derived a way to use the existing code base and 
adapt it to work more like an actual A'-body simulation, where the 
effect of scatterings on the BBH orbit is taken into account, giving 
the orbital evolution. 

In this paper we present a new scheme we developed to sim- 
ulate the inspiral of the BBH from the galactic scale to the hard 
binary scale. Our method is based on imposing conservation of to- 
tal energy and angular momentum, instead of directly summing 
the forces of individual stars on the BHs, and lies between scat- 
tering experiments and A'-body simulations. Studying this type of 
systems with scattering experiments is significantly faster than a 
full A'-body treatment and is appropriate when the bulge has re- 
laxed. The phase space stellar distribution is followed accurately 
and a more realistic number of stars can be included. Thus, we are 
able to run many simulations and probe a large range of parameters 
under reasonable physical assumptions, as well as the convergence 
of the solution, and compare the inspiral timescale, the stalling ra- 
dius, and the eccentricity evolution to earlier calculations. Since in 
our method we obtain the forces on the BHs, it is also possible 
to test analytical expressions for dynamical friction. The original 
treatment by Chandrasekhar (1943) is for a homogeneous back- 
ground, therefore we compare our results to Just et al. (201 1) who 
investigate dynamical friction in power law cusps. Our solution ex- 
tends beyond the range of validity of their formula and into the hard 
binary regime of the BHs. 

We studied the evolution of both an equal mass binary and a 
10:1 mass ratio binary under different initial conditions. We found 
that both cases presented stalling of the semi-major axis, but in the 
unequal mass case, eccentricity tended to grow towards unity on 
timescales well below the Hubble time. 

In Section 2 we give a general mathematical formulation of our 
simulation scheme. The model for the specific simulations we per- 
formed in this work is described in Section 3 while a technical 
description of the scheme appears in Section 4. In Section 5 we 
present the results of all BBH simulations and discuss their phys- 
ical implication. In Section 6 we discuss dynamical friction, com- 
pare to earlier results, and show that our code is applicable in very 
large radii, thus potentially helpful in future studies of phenomena 
related to dynamical friction. Finally, we give a short summary in 
Section 7. 



2 MATHEMATICAL FORMULATION 
2.1 Motivation 

The two basic ideas of our solution are the separation of timescales 
and the balance of energy and angular momentum between the BHs 
and the stars. There are three timescales in the system correspond- 
ing to changes in: 

(i) stellar orbits following close encounters 

(ii) the BHs' orbits 

(iii) the background stellar potential 

In more detail: a star significantly changes its original orbit dur- 
ing a close encounter with a BH, but the BH's path is only slightly 
perturbed as the force exerted on it is due to many small 'scatter- 
ing' events; the change of the background stellar potential is due to 
the collective response of stars to the perturbation, which evolves 



on the dynamical timescale. A yet longer timescale would be of 2- 
body relaxation, which is longer by a factor of ~ O.lN/lnN (Bin- 
ney & Tremaine 2008) and is expected to be well above the Hubble 
time for a real galaxy core. However, relaxation time is potentially 
significantly shorter if the dominant relaxation mode is not 2-body 
relaxation. 



The most basic kind of an A'-body simulation uses very small 
steps, after each the vector forces exerted by all the field stars (par- 
ticles) are summed up and applied to the BHs in the next step; each 
particle is propagated the same way. Energy and angular momen- 
tum are globally conserved (within a given numerical error toler- 
ance) since it is a closed system. If N is not large enough, noisy po- 
tential and unrealistically massive stars lead to spurious relaxation. 
If the time step is too large, close encounters cannot be resolved. 



Our scheme uses a 'large' At, which we call an interval, within 
which there is one or more actual steps of the ODE solver. Let BH 
number i (where i £ { 1 , 2}) propagate one interval between times 
fo and t\ = fo +Af from some vector position S;(*o) to s,-(fi); we use 
the term segment for the path length. After each interval, the ener- 
gies and angular momenta of all the particles are summed up. If no 
work was done and no torque exerted due to the background stellar 
potential, then the energy AE* (angular momentum AL*) gained by 
the stars during some interval must be equal and opposite to the en- 
ergy AEbh (angular momentum ALbh) lost to the BHs during this 
interval. For simplicity, we use a static and spherically symmetric 
model for the stellar potential, so these requirements are automat- 
ically fulfilled. The scheme can accommodate more complicated 
models as well: the demand for spherical (or more generally, ax- 
ial) symmetry, which is necessary if one assumes that change of a 
star's angular momentum is only due to interaction with the BHs, 
can be relaxed if the torque component due to the stellar bulge is 
considered separately; a slowly varying potential can also be taken 
into account if the proper adjustments are made to the code. 



The evolution of a BBH is dominated by close and fast encoun- 
ters with the field stars and the evolution of the stellar gravity field 
is only a secular effect. Thus, the part of AE+ due only to the change 
in stellar potential is negligible after each interval, justifying the 
use of a static model. However, at late times, the initial stellar po- 
tential is no longer consistent with the actual spatial distribution. 
The static potential assumption is still reasonable as far as the BHs 
are concerned, since by the time that any significant evolution of 
the background potential has taken place, the BHs will have fallen 
deep enough in the potential well, where the dominant force is the 
other BH's gravity rather than the background stellar potential (so 
that the exact shape of the potential well does not matter). Stars 
further out orbit in a 'wrong' potential, but as long as approximate 
spherical symmetry is preserved, their interaction rate with the BHs 
should not be significantly influenced by this (and during a close 
encounter, the background stellar potential is of course unimpor- 
tant). 



Assuming additionally that the interval is short enough so that the 
forces on the BHs due to the background stars do not vary signif- 
icantly, we perform scattering experiments in each interval to find 
AE* and AL*. Using the simple algebra described below, we find 
the average forces on the BHs over this interval. We also assume a 
purely planar motion of the BHs. Thus, unless stated otherwise, by 
'angular momentum' and 'torque' we mean only the z component 
of these vectors. 
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2.2 Conservation laws 



The basic conservation equations are 



A£ B h = -A£* ■■ 



si(fi) 

Sl(fo) 



Fidsi 



Sl('l 



s 2 ('o) 



F 2 -ds 2 , 



(1) 
(2) 



AL BH = -AL + = / Tidf + / T 2 dr, 

J/ ./fo 

where F,- and T; = (r,- x F,-) • z are the force and torque, respectively, 
exerted on BH number i by the stellar population. Let us write F, 
in the following form: 



*i = -fi*i+fiti, 



(3) 



where v; and u, are the unit vectors parallel and perpendicular to 
the velocity of the BH. Note that / is 'drag like' and directed op- 
posite to the velocity vector. If an object moves through a uniform 
background, symmetry dictates that the mean force would be an- 
tiparallel to the velocity vector. However in realistic environments 
there must be also a force due to the inhomogeneities of the back- 
ground, so a perpendicular force component is required. In Carte- 
sian coordinates: 



— x+^y, 

Vj V; 

— * y 

V; V; 



(4) 
(5) 



To simplify equation (1) we write ds,- = \,dsj and thus 

s,('i) rutin) . 

F, ds, = / (-/,v, +/,u,) -v/dj,- 

Sf('o) • / S/(fo) 



/•s,(«i) 

fi / ds; = -fiSi. 



(6) 



The perpendicular force component / disappears due to the dot 
product, while / is assumed constant along the path and can be 
taken out of the integral, which defines S (which is simply the 
path's length). Finally, equation (1) for the energy becomes 

AE i , = f 1 S 1 +f 2 S 2 . (7) 

To simplify equation (2) we write Ti in Cartesian coordinates: 

fi fi 

T; = — - {xiVty - yiVix) - — (xiVjx + y,-Vfy), (8) 

Vi Vi 



and thus 



f ' Tidf = -fi [ 



dt-ft I -At 



to v i 



= -m-fiQi- 



(9) 



The above integrals define V and Q. Equation (2) for the angular 
momentum becomes 



al* = fiVi +/i Qi + fiPi + hQi- 



(10) 



If we evolve the BHs between to and t\ under their mutual gravity 
alone, energy (and angular momentum) will be conserved along the 
produced orbital segment: AEbh = 0- Stars scattered on this orbital 
segment will not their conserve energy: AE* ^ 0; so energy is also 
globally not conserved between to and t\ . By solving equations (7) 
and (10) for the (non conserving) forces, the BHs can be evolved 
again in this time interval, under the additional forces, producing 
an orbital segment for which A£bh = ~ &E±- The revised orbital 
segment is only slightly different from the original, since the addi- 
tional forces are much smaller than the forces exerted by the other 
BH and the background potential. Stars scattered on the revised or- 
bit will have a slightly different AE*. This process of altematingly 
evolving the BHs and the scattering of stars can be repeated until 
converges is achieved. 

Since AE* and AL* are directly obtained from the scattering ex- 



periments and the six integrals (calligraphic letters) are calculated 
from the orbital segments, equations (7) and (10) are a linear system 
of two equations with four variables: f\, / 2 , f\ and / 2 . It is worth 
noting that without the perpendicular force, which is expected to be 
negligible in the standard dynamical friction formalism, it is impos- 
sible to conserve energy and angular momentum simultaneously. In 
particular cases, additional constraints give an exact solution as ex- 
plained in the following section. 



2.3 Finding the forces 

2.3.1 Symmetric motion 

In this case, the masses are equal and the initial conditions are sym- 
metric with respect to the centre of the system. The motion of one 
BH mirrors that of its companion, so that the path integrals are 
equal for the two BHs (the index is therefore dropped) and the 
forces acting on them must also be equal due to symmetry. The 
solution is 



f = — A£*, 
' 25 



/ = 



2SQ 



[AL+S-AE+P). 



(11) 
(12) 



There is no solution for Q = (but S ^ is guaranteed by equa- 
tion 6). In points where the BHs' velocity is purely tangential, / 
is parallel to the radius vector and thus exerts no torque (and work 
is never done by /). At these points, / is free but /; is overdeter- 
mined (must change both £bh and Lbh by the specified amounts). 
Segments which are symmetric about such points have Q = 0. 



2.3.2 High mass ratio 

In this case one mass is much larger that the other, and is assumed 
to sit motionless at the centre of the system. The forces / and / 
act on the secondary BH only. The integrals S, V and Q are also 
calculated for the secondary only. The solution for Q ^ is 

1 



/ = 



S 
1 

SQ 



(AL+S-AE+V). 



(13) 
(14) 



2.3.3 General solution 

The results for the limiting cases described above motivate us to 
look for solutions of the form 



h = ^A£„ 

l-a 
h = ^^A£, 



h 
h 



SiQi 

l-a 



(AUSi-AE*Vi), (15) 
(AL*S 2 -A£*P 2 ). (16) 



5 2 Q 2 

There is a mathematical solution for every a, but we know from 
the limiting cases that ^ a < \. The value of a is a function of 
the mass ratio q, but it may also be dependent on other factors, such 
as the local stellar densities at the instantaneous position of either 
BH. In the case of a high mass ratio, one BH is almost stationary; 
therefore a must approach zero faster or at least as fast as <Si . The 
force f\ does not have to approach zero, but the acceleration f\/M\ 
does. 



3 MODEL 
3.1 Units 

For reasons of consistency with our previous work (Meiron & Laor 
2010), we use a unit system in which mass is measured in units 
of the primary BH's mass, velocity is measured in units of 4rr and 
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G = 1 (where a and G are the stellar velocity dispersion and the 
gravitational constant, respectively). The hard binary separation is 
defined as 



a h - 



q GM. 

l+q 4o2"' 



(17) 



where q ^ 1 is the mass ratio of the secondary and primary BHs 
and M. is the mass of the primary. In our units, the hard binary 
separation of an equal mass binary (q = 1) is 2. 

The base units are therefore scalable by M. and o". Only one 
parameter is required if we also use the M-a relation (e.g. Gultekin 
et al. 2009). The units of length, time and velocity units and their 
scaling, using M. and a, are 



[L] = j^GM,g~ 2 = 0.77 Mg' 53 pc, 
[T] = ^GM.e- 3 = 1000 Ml 
[V] = 4(T = 750 Mg' 24 km s" 



0.29 



yr, 



(18) 
(19) 
(20) 



where M% is the physical mass of the primary BH in units of 
10 8 Mq. Note that in all simulations we used a stellar velocity dis- 
persion of 0.25 velocity units. 



3.2 Bulge Properties 

In all our simulations, stars are initially distributed in a singular 
isothermal sphere and follow a Maxwell-Boltzmann distribution 
with ID velocity dispersion a. To avoid the non-physical diver- 
gence of the potential, we assume a core structure: 



P{r) 



Po 
Po 



r> h 



(21) 



where h is an arbitrary break radius set to 1 and po = o" 2 /2nGh 2 . 
The expression for the gravitational potential (or the bulge poten- 
tial) derived from the above density is 



*bul g eM = i ^ 



[f- +21n(£)-l 



r^h 
r> h 



(22) 



Koopmans et al. (2009) found that massive elliptical galaxies, 
within their effective radii, are well approximated by a power law 
ellipsoid with an index of —2. Genzel et al. (2003) also found that 
the density of the nuclear star cluster of the Milky way can be de- 
scribed by a broken power law with index of —2.0 ±0.1 down to 
0.38 pc. This is however not appropriate within the BH sphere of in- 
fluence, if dynamically relaxed, where the equilibrium distribution 
(radius independent mass and energy flow) is the famous Bahcall 
& Wolf (1976) cusp of p « r~ 7 / 4 (where only one mass species 
is present, cf. Bahcall & Wolf 1977, Hopman & Alexander 2006). 
Unrelaxed clusters around adiabatically growing BHs are expected 
to have shallower slopes (Young 1980) or steeper slopes in the case 
of rotating systems and non-isothermal clusters (Lee & Goodman 
1989). 

In the case of a minor merger (equivalent to q = 0.1), the struc- 
ture of the more massive galaxy does not change significantly and 
equation (2 1 ) likely represents the stellar environment seen by the 
secondary BH after its parent galaxy is absorbed. This picture is 
somewhat naive in the case of a major merger (q = 1), but violent 
relaxation (Lynden-Bell 1967) due to the rapidly varying poten- 
tial in a newly merged galaxy causes widening of the stellar en- 
ergy distribution and is analogous to relaxation by collisions in a 
gas. This process tends to drive galaxies towards a universal steady 
state (Syer & White 1998). A nice demonstration of this appears in 
the N-body merger simulations of Milosavljevic & Merritt (2001), 
who find an r~ 2 density profile at the time the equal mass binary 
becomes hard, which extends down to the scale of the binary sepa- 
ration. 



Throughout this work we assume that the BHs are 'naked', that 
is, do not carry clusters of bound stars. In the equal mass case, 
the stellar mass bound to a single BH once the binary becomes 
hard does not exceed 10 per cent of its mass (assuming that the 
cluster also has a power law density profile with index of -2 and 
normalization based on the M-o relation). This will most likely 
not affect the late time evolution, but will surely affect the early 
inspiral phase, and the exact inspiral time, but will probably not 
have a major effect. 



4 ALGORITHM 
4.1 Equal masses 

Below we provide a technical description of the application of the 
technique described in Section 2 for an equal mass binary simula- 
tion. The slightly different procedure for the 10:1 mass ratio case is 
discussed in the following Section. 

First, a realization of a singular isothermal sphere is produced 
up to a cutoff radius of R max , with two equal mass BHs placed 
on the x-axis at x = ±/?oi at this initial distance the BHs are still 
unbound to each other and their inspiraling orbits are governed by 
the bulge's gravity and dynamical friction. The simulation duration 
is divided into (equal) intervals At, after each the stellar force acting 
on the BHs is updated. The actual time steps of the ODE solver are 
smaller if necessary, so close encounters can be resolved. Within 
each interval we do the following: 

(i) Symmetrically advance 1 the two BHs from t t to under 
their mutual gravity, and the effective forces exerted by the stars. 
The calculated antiparallel / and perpendicular / embody the stars' 
pull on the BHs, so the bulge potential (equation 22) should not be 
considered additionally. This gives a short orbital segment which is 
stored in memory. 

(ii) Advance each star from f; to under the influence of the 
two BHs and the bulge potential. The motion of the BHs in this 
interval is already set from the previous stage, so this is essentially 
a short scattering experiment. This is the most computationally de- 
manding stage of each iteration, but easily parallelized. 

(iii) Sum up the energies and angular momenta of all stars and 
subtract the values from the previous iteration to obtain Afi* and 
AL*. 

(iv) Use the BHs' path to obtain S, V and Q. 

(v) Calculate / and / from equations (11) and (12), to be used 
in the next At interval. 

Unless otherwise indicated, our simulations end at t = 10000 
time units (or 10 7 years for a 10 s Mq primary); this is equivalent 
to ~ 1 000 revolutions after the binary becomes hard in most equal 
mass simulations. Integration of a single star is terminated prema- 
turely in three cases: if it reaches a distance of Rmax + 10 length 
units from the centre of the system, reaches r^nl = 10~ 3 from ei- 
ther BH or takes more than 35 integration time steps to complete 
the interval of At = 0.1 (see Section 4.3). In the first case the star 
is considered to have escaped the system (or diverged); the extra 
10 length units beyond R mdx are an arbitrary 'padding' required for 
technical reasons. The second case represents tidal disruption of the 
star (the orbit is then said to have crashed). The true tidal disruption 
radius is up to two orders of magnitude smaller than our Ptidali but 
this choice has a negligible effect on the BHs because the rate of 
crashing stars is negligible compared to the rate of diverging stars. 
The choice of Z? ma x is rather arbitrary, and is chosen to be large 



1 The word 'advance' in this context means solve the equation of motion by 
means of Runge-Kutta method of order five with adaptive step size control 
(e.g. Press et al. 1992). 
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enough to minimize its effect on the results, as discussed in Section 
5.2. 

Notably, the forces on the BHs are one interval retarded. Thus, 
the interval duration must be short enough in order not to break 
conservation of energy and angular momentum. In all simulations 
we chose to use equal intervals of At = 0.1, and since simulations 
done with this choice both ran reasonably quickly and performed 
well in terms of conservation, we did not thoroughly investigate 
changing At. Certainly making the interval length longer or adap- 
tive can significantly speed up the simulations. 

The code outlined above is very simple in terms of decision mak- 
ing; the bulk of CPU effort is made to individually advance stars 
(scattering experiments) with only a small non-parallelizable over- 
head between intervals. A typical simulation with a large realiza- 
tion (N = 5 x 10 6 ) would run for ~ few days on a medium strength 
personal desktop computer. 

4.2 Mass ratio 

As noted in Section 2.3, the force acting on the BHs can be uniquely 
found from the energy and angular momentum differences only in 
the cases of equal masses and high mass ratios (where it is possi- 
ble to assume that the primary BH is fixed at the centre). However, 
when the velocity is perpendicular to the radius vector (i.e. at peri- 
centre or apocentre), the values of Qj approach zero, and thus /; 
are left out of the coupled equations (10) and (7), and cannot be 
solved for. At the same time, /; are overdetermined as they have to 
compensate for both energy and angular momentum changes. Even 
in intervals that contain an apsis, there is usually a solution unless 
the orbital segment happens to be symmetric about the apsis. Since 
only a few segments are affected, the effect on the BBH orbits is 
not significant. 

In the case of equal masses, no special treatment of these seg- 
ments was necessary; but in the case of very high mass ratios (i.e. 
q <C 1) the finite numerical fluctuations are larger due to the fact 
that the secondary BH is less massive. In these cases the total an- 
gular momentum would discontinuously drop at an apsis, and the 
otherwise smooth path would suddenly break at this segment. We 
compromised on simulating a 10:1 binary, and utilized two tech- 
niques to improve accuracy. First, each step was performed twice: 
after calculating the frictional force, the stars and BH were reset 
to their original positions, and advanced again with the newly cal- 
culated force on the BH. Second, we attempted to compensate for 
the accumulated error: instead of calculating the force components 
using purely AE* and AL* of the last At interval, we added their re- 
spective accumulated errors (the correction term had weight of 0. 1 
per cent). Those adjustments dramatically improved the accuracy of 
the 10:1 simulations, but unfortunately we did not yet overcome all 
the technical problems associated with reliably simulating a higher 
mass ratio inspiral from ~ 50 pc down to stalling separation. 

4.3 Quality Assurance 

If the routine that advances the stars has a bug or is just not accurate 
enough, the values of and AL* obtained in each interval will 
be faulty; but since the frictional force is calculated in such a way 
that would compensate for any change of E and L in the stars, the 
bug might remain undetected. Thus, an important validity test of E 
and L conservation in the code is inherently not available. 

We can, however, test the ODE solver for a similar problem, and 
infer that our stellar orbits are at least well behaved in the original 
problem. Assume a BBH with a circular orbit of constant radius R; 
the BBH orbit is now decoupled from the stars and each stellar orbit 
can be integrated separately. This is the restricted 3-body problem 
plus a spherically symmetric potential; as in Meiron & Laor (2010), 
the energy in the rest frame of the BHs (the frame which rotates 



with the same angular frequency) is the only conserved quantity. 
The Jacobi integral is a constant of motion related to the rest frame 
energy by Cj = — IE. 

We performed a number of tests where the BBH was in a cir- 
cular orbit with a radius in the range of 1 $J R ^ 40 length units. 
The stellar model was the same as described in Section 3.2 and the 
orbits were evolved for t = 10000 time units. It was found that Cj 
is well conserved for the great majority of stars: only 0.2 to 1.6 per 
cent of orbits were cast off as rogue orbits, exceeding 35 steps per 
At = 0.1 interval; the rest had an average |AC/| of (4 to ll)xl0 -5 
energy units (Cj is typically of order unity). When applying a short 
softening length of 0.04 (corresponding to Quinlan & Hernquist 
1997), 0.2 to 1.7 per cent of orbits were eliminated by the same cri- 
terion and the rest had an average \ACj\ of (5 to 17) x 10~ 5 energy 
units. Thus, we did not apply softening in the actual simulations. 
This type of number of steps filter against rogue orbits was found 
to work better than putting a lower limit on the allowed step size. 

Accuracy may be improved simply by integrating close encoun- 
ters with smaller error tolerance, or by employing more elaborate 
techniques such as regularization. The small number of stars which 
are lost due to the ODE solver's inability to handle them does not 
have a significant systematic effect on the BBH orbital evolution. 



4.4 Conservation of E and L 

As noted previously, total energy and angular momentum conserva- 
tion is not a built-in requirement of the method, but rather indicates 
a successful transfer of energy and angular from the BBH to the 
stellar population. In the upper panels of Figure 1 we show energy 
of the BBH (solid black line), the stellar population (solid red line) 
and their sum (dashed blue line); the lower panels is the same but 
for the angular momenta. The left panel is for an equal mass simu- 
lation while the right panel is for a 10:1 mass ratio simulation. 

The quality of energy conservation is attested by the absolute 
difference of total energy between the beginning and the end of the 
simulation; it is also useful to look at the fluctuations in total en- 
ergy, a very crude estimate of which is the amplitude of the largest 
fluctuation (this is not a very good measure since fluctuation am- 
plitude can be high in some parts of the orbit mild in others). These 
two quantities have dimensions of energy, and it is most sensible to 
normalize them with respect to the absolute difference in the energy 
of stars (or the BHs). Thus, for a simulation ending at t — T: 

sf(T) = [£ tot (r)-£ tot (0)]/|^(r)-£ + (0)|, (23) 
ef (r) = max(|£ tot (f) - (E tot >|)/|£*(r) -£*(0)|. (24) 

We similarly define ef and ef for the angular momentum. All these 
e parameters must be very small. If the quantity that is supposed to 
be conserved has some trend, then usually £j = £2, otherwise it only 
fluctuates around its mean value and £j < £2- 

For an equal mass simulation with N = 5 x 10 (discussed in de- 
tail in Section 5.1) case at T = 43200, we get ef ,e[ < 10~ 6 and 
e 2 = e 2 = ^ x 10 4 - F° r a 10: 1 mass ratio simulation (discussed in 
detail in Section 5.4) at T = 21 800, we get ef = 4 x 10~ 6 and ef = 
6 x 10~ 5 , with a decreasing trend in total energy: ef = ef = 0.008. 
The trend in total energy begins at t ~ 14000 and is probably due 
to the fact that there is an accumulated inaccuracy in the solution 
of the BHs' equations when they are so close together at pericen- 
tre. At t = 14000 we get ef = 8 x 10~ 5 and ef = 4 x 10~ 4 with 
values for the angular momentum similar to the end of the simu- 
lation. Thus, the numerical effect previously discussed is unrelated 
the eccentricity growth observed in this simulation. 



6 Y. Meiron and A. Laor 



13.2 I 1 1 1 1 1 I ~| 11.5 




Figure 1. Conservation of energy and angular momentum in the case of equal masses on a circular orbit (left panel; same simulation as shown in Fig. 4), and 
the case of 10:1 mass ratio with initial eccentricity of 0.2 (right panel; same simulation as shown in Fig. 8). The dashed blue lines are the system's total of 
the conserved quantities per unit mass; the solid red line represents the stellar population's part and the solid black line represents the BHs' part. The BH-star 
interaction term is in the stellar part, so the BHs' energy in this figure is only the sum of their kinetic energies and mutual potential energy. 
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Figure 2. A representative sample of simulation results showing the radius 
as a function of time. The red and blue lines are for an initially circular and 
eccentric equal mass binaries, respectively; magenta and cyan are the same 
for a mass ratio of 10:1; the green line is the same as the red line but with 
a larger cutoff of the stellar bulge. The thickness of the line is indicative of 
eccentricity. While the semi-major axis stalls in all simulations, eccentricity 
does not reach a steady state in the unequal mass cases. 



5 RESULTS 

We performed a total of 58 simulations of an equal mass binary 
and a 10:1 mass ratio binary, varying N, Rq, R m ax and the initial 
eccentricity eo. In Figure 2 we shows r{f), for selected simulations. 
In the 10:1 simulations, the primary BH is fixed at the centre, so r 
is the BBH separation; in the equal mass simulations r is half the 
separation. As can be seen in the figure, the semi-major axis stalls 
in all simulations but eccentricity (indicated by the thickness of the 
lines) does not reach a steady state in the unequal mass cases. These 
results are discussed in more detail below. 



5. 1 Number of Particles 

Here we show two things: how the results converge with increasing 
number of stars N, and how the results depend on the specific re- 
alization of the stellar distribution. We present seven pairs of sim- 
ulations with N between 50000 and 350000; for each N the two 



simulations have a different random seed, so that the stars have dif- 
ferent initial positions and velocities, but are drawn from the same 
distribution. All simulations are of an equal mass binary starting 
at Ro = 60 with the local circular velocity; the cutoff radius of the 
stellar sphere is ^ max = 200. 

In Figure 3 we show the 'final state' (i.e. at t = 10000) semi- 
major axes a (equivalent to separation, in the equal mass case), and 
eccentricities e as functions of N; these two numbers are the best 
way to appreciate differences between similar simulations. From 
this small sample, it is apparent that the effect of changing N in 
this range is comparable in magnitude to that of changing the re- 
alization. In this set of simulations, the range of semi-major axis 
values is 0.07 length units or ~ 5 per cent of the sample's average 
a; the eccentricities are small and in the range 0.03 < e < 0. 1. Thus, 
increasing N beyond 100000 (within R max = 200) is unnecessary 
for this level of accuracy and following tests are made using this 
number. 

The average of a in these simulations is 1.43 (in physical units, 
for 10 8 Mq BHs, this is equivalent to 1.1 pc); this is approximately 
30 per cent below the hard binary separation of 2. Merritt (2006) 
suggested the following formula for the stalling separation: 

"stall 



0.2 



(l+<7) 2 



(25) 



where r' h is the radius containing a mass in stars equal to twice the 
combined mass, or M(rl) = 2A/,(1 +<?), at the time of stalling. 
We measured the accumulated mass in the above simulations and 
got that in all of them r, = 37.5 with a very small spread. By sub- 
stituting this information and q = 1 into equation (25), one gets 
fl stall = 1-88. It is important to note that in the Merritt (2006) sim- 
ulation there was no actual stalling of the BBH, and a sta u was es- 
timated as the value of a in which a clear change in the hardening 
rate took place; r 1 , was determined at the time when this change 
occurred. In our simulations both r' h and a sta u were determined at 
t — 10000, which is well after the hardening rate has dropped. 
Also, equation (25) was calibrated by Merritt using simulations 
with mass ratios 0.025 q 0.5. 

To test the convergence of the results, we also performed a sin- 
gle simulation with significantly more particles (N = 5 X 10 6 ) and 
longer duration (f = 43200); we show the inverse semi-major axis 
as a function of time in Figure 4. As seen in the figure, there is 
still some slow evolution of the semi-major after / = 10000. At 
t = 40000 the value of a is 1.35, which is 5.6 per cent lower than 
the value at t = 10000 and only one per cent lower than the value at 
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Figure 3. The semi-major axes a and eccentricities e at t = 10000 as func- 
tions of the number of particles N, for equal mass and initially circular bi- 
naries. For each N, the two circles represent two different realizations of 
the same stellar population. The effect of changing N in the range tested in 
this set is comparable in magnitude to that of changing the realization, and 
no trend is seen. The simple averages for the entire set are a = 1 .43 ± 0.02, 
corresponding to a stalling separation of ~ 1.1 pc for two 10 s Mq BHs, and 
e = 0.06 ± 0.02 which indicates nearly circular orbits. An additional run 
with N = 5 x 10 s is not shown here (see text). 
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Figure 4. The inverse semi-major axis of an equal mass binary simulation 
that was performed with N = 5 x 10 6 particles (compared with N ^ 3.5 x 
10 5 for the simulations in Fig. 3) and for a much longer duration. The value 
of a at t = 40000 is 1.35, which is 5.6 per cent lower that the value at 
t = 10000 and only one per cent lower that the value at t = 30000. Cf. 
figure 1 of Berczik et al. (2006) 

t = 30000; the decay rate at the end of this simulation is d = 10~ 6 
velocity units, equivalent to ~ 10~ 9 pc yr _1 for a 10 s Mq primary. 
In this specific run, eccentricity was especially low at e < 0.01. 

The orbits produced in the other simulations in this set (with 
N ^ 3.5 x 10 5 ) are very similar to the orbit shown in Fig. 4, re- 
gardless of the number of particle (cf. Berczik et al. 2006, spherical 
and triaxial cases). Note that a large N is required only to minimize 
the statistical fluctuations, our scheme is not subject to an artifi- 
cial stellar relaxation mechanism, which requires a very large N to 
overcome. 



5.2 R max and R 

Since we only simulate the spherical component and not a full 
merger, the initial and boundary conditions need be assumed: the 



Figure 5. The semi-major axes a at t = 10000 as a function of the cutoff 
radius R mdx , for equal mass and initially circular binaries. For each of S max , 
the tree data points represent BHs launched from different radii Ro- The in- 
crease of a with Ro seems to saturate, while the dependence on i? max is very 
weak in the tested range. Eccentricity in this simulation was not correlated 
with either of the tested parameters, its simple average for the entire set is 
e = 0.06 ±0.03 which indicates nearly circular orbits. 

initial distance of the BHs from the centre, Rq, and the cutoff radius 
of the stellar sphere, R mdx . Here we present six trios of simulations 
with R ma x between 120 and 520, each three simulations are of an 
equal mass binary starting at Rq = 30, 60 and 90 with the local 
circular velocity. The number of particles is N = 500 x R mSK , so 
as to keep the particle density profiles of equal normalizations (but 
different cutoff radii) in all simulations; this is due to the fact that 
the number of particles (or mass) grows linearly with radius in an 
isothermal sphere. 

In Figure 5 we show the final state semi-major axes a (stalling 
separation) as a function of R max . From the results discussed in 
Section 5.1, a characteristic error of ~ 0.04 units on the semi- 
major axis can be attributed to a specific realization and number- 
of-particles statistical fluctuation. The simulations with Rq = 30 
give systematically lower values for a; the values for Rq = 60 do 
appear to be systematically lower than the Rq = 90 simulations, but 
the two sets are within the errors of each other. There appears to be 
a trend of decreasing values of a with increasing A ma x up to ~ 300, 
which is to be expected both because there is a larger supply of par- 
ticles that can interact with the BHs, and because the potential well 
is deeper and stars that have already interacted with the BHs have a 
larger probability to fall back to the centre and interact again. Nev- 
ertheless, we see from this small sample that these effects are weak 
and comparable in magnitude to those discussed in the Section 5.1. 
In this set of simulations, the eccentricities are also small and in the 
range 0.02 < e < 0. 11; no correlation of eccentricity was observed 
with either Rq or ^ ma x- 

5.3 Mass Ratios 

As noted in section 4.2, the simulation of an unequal mass binary 
is somewhat different in nature. Thus, when studying the evolution 
of a 10:1 binary, we followed the analysis of Section 5.1 and per- 
formed a number of different tests with increasing number of par- 
ticles and different realizations. Here we present four pairs of sim- 
ulations with N between 100000 and 400000; for each N the two 
simulations have a different random seed, so that the stars have dif- 
ferent initial positions and velocities, but are drawn from the same 
distribution. All simulations are of a binary with a 10:1 mass ratio 
starting at Rq = 60 with the local circular velocity; the cutoff radius 
of the stellar sphere is R mlix = 200. 

In Figure 6 we show the final state semi-major axes a and eccen- 
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Figure 6. Same as Fig. 3 but for a 10:1 binary. The semi-major axes are 
consistent within this sample, the simple average is a = 0.241 ± 0.005, cor- 
responding to a stalling separation of ~ 0. 19 pc for for a 10 7 M s secondary 
around a 10 s Mq primary. Eccentricity varies greatly in this sample, see 
text. 



Table 1. A list of all 10: 1 binary simulations. The columns from left to right 
are: the number of particles, initial eccentricity, eccentricity at t = 10000 
and the estimated timescale to reach e = 1 (see text for definition). 
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tricities e as functions of N. As in Section 5.1, changing N in this 
range produces no apparent effect on the stalling radius. Within this 
set of simulations, the semi-major axes are consistent; the sample's 
average a is 0.241 ± 0.005. We measured the accumulated mass 
in the above simulations and got that in this case r' h = 16.8 with 
a very small spread. By substituting this information and q = 0.1 
into equation (25), one gets a sta u = 0.28. The caveats of using this 
equation were explained in Section 5.1. If we recalibrate equation 
(25) using our two values of q, the prefactor is lowered from 0.2 
to approximately 0.16; it is even somewhat smaller considering the 
fact that the true value of a sta u can be ~ 6 per cent lower than its 
measured value at t = 10000 (see Figure 4). 

In contrast with the equal mass case, the eccentricities do not 
reach a steady state value. In one of the N = 10 5 simulations, the 
eccentricity increases very rapidly after the binary becomes hard. 
At / = 9400 the pericentre distance reached 10~ 3 and the simu- 
lation is terminated. In the rest of the simulations, The eccentric- 
ity values range between 0.05 and 0.28 at t = 10000, however, in 
all but one of these simulations eccentricity is slowly increasing. 
We calculate T e ^\, a very rough estimation for the time of grav- 
itational wave regime, by fitting e(t) with a linear function in the 
range 6000 < t < 10000 and continue it to e = 1. At this time span 
the rise in eccentricity is approximately linear, however this trend 
breaks at approximately e = 0.85, so in fact T e ^\ can underesti- 
mate the time to the gravitational wave regime by some ~ 20 per 
cent. In Table 1 we show the results for T e ^\ for this simulation set 
and also for the initially eccentric runs. The values range between 
approximately 4 x 10 4 and 3 x 10 5 time units, or on the order of 
10 8 years for a 10 7 M© secondary around a 10 8 Mq primary. 



5.4 Eccentricity growth 

Here we tested how the initial eccentricity affects the results. We 
performed simulations with four different initial eccentricities be- 
tween 0.1 and 0.4 with semi-major axes as in the circular simula- 
tions of Sections 5.1 and 5.3. The initial eccentricity eo is that of 
a BBH with the same initial conditions and that moves in the ini- 
tial potential (equation 22) with no friction; the orbit is not really 
an ellipse since the potential in not Keplerian, so eo corresponds 
to the mean orbital eccentricity, defined as the difference between 
the maximum and minimum separations divided by the major axis. 
For each eo, there are two simulations (with two different realiza- 
tions) for an equal mass binary and two for a 10:1 binary. The 



number of particles is N = 2 x 10 5 with cutoff radius i? ma x = 200; 
we performed a single simulation with significantly more particles 
(N = 5 x 10 6 ) for the e = 0.2 case. 

In Figure 7 we show the final state semi-major axes a and ec- 
centricities e as functions of eo. This figure also includes four data 
points with eo = that have already been presented in Figs. 3 & 
6. The stalling separation is independent of initial eccentricity in 
the tested range, the sample's average a is 1.42 ±0.04 for the equal 
mass simulations and 0.246 ± 0.010 for the unequal mass simu- 
lations. As in the other tests we performed for equal mass bina- 
ries, the final orbits are very much circular; the eccentricities are 
in the range 0.01 < e < 0.1 despite the initially significant eccen- 
tricity. However, the final eccentricities in the 10:1 cases do ap- 
pear to be generally correlated with eo- As in Section 5.3, here too, 
eccentricities are still increasing when the simulations terminate 
at / = 10000; the eccentricity timescales T e _>\ for the 10:1 sim- 
ulations are shown in Table 1, which shows that the eccentricity 
growth rate is related to eo. 

Sesana (2010) found that eccentricity growth is generally mild 
for equal mass binaries with very small initial eccentricities, but 
also that initial eccentricity eo > 0.3 leads to very high peak eccen- 
tricity almost regardless of the other system parameters. Although 
all our equal mass simulations end up in nearly circular orbit, this 
is not inconsistent with Sesana: while in the latter work the BHs 
are launched from within their radius of influence, in our simu- 
lations the BHs are launched from much further out, in the dy- 
namical friction regime. In our simulations, an equal mass binary 
becomes less eccentric as it inspirals from r > 60. For one of the 
eo = 0.4 simulations, the rapid eccentricity decrease ceased at about 
r = 1.4 (approximately twice the stalling radius) where the value 
was e = 0.035. At r = 4 (equivalent to the initial radius of Sesana's 
equal mass binaries) the eccentricity was 0.12; no significant ec- 
centricity growth occurred in Sesana's equal mass simulation with 
initial eccentricity of 0.1. While it has been shown by Dotti, Colpi 
& Haardt (2006) that BBHs lose memory of their initial eccentric- 
ity if they corotate with a massive gaseous disc, studies of eccentric 
orbits of hard binaries is motivated for the purely stellar dynamical 
case by the theory of linear response for dynamical friction (Colpi, 
Mayer & Governato 1999). However, this theory is derived from a 
first order perturbative expansion and is not applicable when close 
encounters dominate the evolution, and the system is not well de- 
scribed by an analytical approximation. 

It is important to note another major difference between our work 
and Sesana (2010) which greatly affects the evolution of the binary 



A conservation-based method for simulating BBHs 9 



1.6 
1.5 
1.4 

a T 
0.3 - 
0.2 

1.0 ~ 

0.8 - 
0.6 - 
0.4 - 
0.2 - 
0.0 



mm m ^ 

H ± ^ 



0.0 



0.1 



0.2 



0.3 



0.4 



Figure 7. The semi-major axes a and eccentricities e at t = 10000 as func- 
tions of initial eccentricity eg, for equal mass (blue circles) and 10: 1 binaries 
(red squares). For each eg, the two data points of each kind represent two 
different realizations of the same stellar population (all simulations have 
discrete values of eg, but some overlapping data points were moved slightly 
to the left and right for graphical reasons). In the 10:1 case, the final ec- 
centricities increase with eg while the semi-major axes are unaffected. The 
equal mass binaries seem to circularize independently of their initial ec- 
centricities, as all results are very similar in both a and e. The red stars 
represent a single 10:1 simulation with N = 5 x 10 6 ; all others shown here 
have A? = 2 x 10 s . 



orbital parameters, mostly the semi-major axis: while in our work 
the loss cone empties, Sesana implicitly assumes that the loss cone 
is always full at r > r; n f. This leads to a very rapid decay of the 
binary separation and quick coalescence due to gravitational wave 
emission. Similarly, Antonini & Merritt (2012) calculated the ec- 
centricity evolution in the case of a very small secondary BH that 
does not affect the stellar density profile. There, eccentricity grows 
because the orbit passes in and out of a flat core, where the star 
are fast and the drag force is much less efficient at pericentre than 
at apocentre. In our simulations, however, the much more massive 
secondary forms a cavity slightly larger than its apocentre, and ec- 
centricity grows where the density is essentially zero. 

In Figure 8 we show the evolution of the semi-major axis for all 
the runs with eo = 0.2 and 10:1 mass ratio. The dotted green lines 
represent the two realizations with 2 x 10 5 particles while the solid 
blue line represents the larger N = 5x 10 6 realization. The first two 
simulations are arbitrarily terminated at t — 10000 while the latter 
is stopped only when the eccentricity reaches 0.99. If we scale to 
physical units for a primary BH of 10 s M@, the end of the simula- 
tion is 22 Myr from its beginning. For this mass scaling, using the 
Peters (1964) formula for orbital decay due to gravitational waves, 
the timescale for coalescence at the end of the simulation is just 1 
Myr. 

The rapid growth of the eccentricity while the semi-major axis 
remains fairly constant indicates a high value for the the dimen- 
sionless eccentricity growth, defined as: 

K= % . (26) 
dln(l/a) 

For single scattering of unbound stars from a fixed background, 
Quinlan (1996) derived a maximal value of K of about 0.3 for mass 
ratio of 16: 1, consistent with previous scattering experiments (Roos 
1981; Mikkola & Valtonen 1992). The value of K for the simula- 
tion presented in Fig. 8 is at least an order of magnitude larger (we 
only roughly estimated the value from the results). The difference 
is probably due to the very different nature of the orbits in the re- 
stricted 3-body problem versus the more realistic model used here. 



We refer an in-depth study of the physical mechanism behind the 
eccentricity growth to a subsequent work. 

It is interesting to compare our results to those of Iwasawa et al. 
(2011), who also got 'runaway' growth of eccentricity while the 
semi-major axis stalled, but one must note the critical differences 
between the two studies. Most notably, Iwasawa et al. used a mass 
ratio of 100:1 while we took only 10:1. Additionally, they used 
an initially very shallow central density profile, p °= r~ 3 / 4 while 
our bulge model (see Section 3.2) was an isothermal sphere, p °= 
r~ 2 . More importantly, the mass of their entire stellar population 
was just 1.25 x 10 9 M ( o, which is 8 times less than the primary 
BH's mass and only 12.5 times more than the secondary's mass. By 
comparison, the total stellar mass in our model is 25 times the mass 
of the primary and 250 times the mass of the secondary. Thus, their 
entire simulation is deep within the primary's sphere of influence, 
where its gravity dominates, while our simulations started with the 
secondary well outside the primary's radius of influence. 

Let us scale our work to theirs by setting Mg = 100 in the scal- 
ing equations of Section 3.1. Iwasawa et al. start their A32k sim- 
ulations (with N = 32768) at Rq = 20 pc within which there are 
less than 2500 particles, their BH stalls at a = 3.9 pc. We initiate 
the secondary BH at r = 637 pc and get stalling at 2.1 pc; there 
are initially 40000 particles enclosed within r = 20 pc in our large 
N — 5 x 10 6 simulation (marked with a star in Fig. 7). The sec- 
ondary BH is 20000 times more massive than a field star in our 
simulation, versus 2600 in the Iwasawa et al. simulation. Their stel- 
lar bulge model is very small compared to ours, with 90 per cent of 
their mass is within Ro = 190 pc, while our model is truncated at 

tfmax= 1800 pc. 

The semi-major axes in the two studies evolve at very different 
rates: while the Iwasawa et al. binary takes 19 million years (Myr) 
to sink from oq = 20 pc to 10 pc in their shallow cusp, our binary 
does the same journey in only 0.18 Myr. However, when applying 
equation (6) of Iwasawa et al. (201 1) (derive from their numerical 
results) for our physical parameters, the timescale for significant 
eccentricity growth is extremely short at 0.78 Myr; the lifetime of 
the system from our simulation is approximately 85 Myr (scaled 
with Mg = 100), which is two orders of magnitude longer. The dif- 
ference might be due to the mass scaling assumed in their formula. 

A final note about precession: in the unequal mass case, the BHs 
exhibit very small precession during the hard phase. In the large 
simulation of Fig. 8, between t\ = 10000 and 12 = 1 1 000 the semi- 
major axis precesses by 0.31 1°. During this period of 1 Myr (scaled 
with Mg = 1), the average semi-major axis is 0.18 pc, and it drops 
by 0.002 pc; the average eccentricity is 0.559 and the growth is by 
0.072. This precession can be produced, for example, by perturb- 
ing the Keplerian potential with a uniform density field of 20,000 
solar masses per parsec cubed. Even with our large number of parti- 
cles, this density corresponds to only four particles enclosed in the 
sphere with radius equals to the apocentre. In the snapshot taken at 
t\ , there was one particle inside this region, and it was probably a 
transient since there are no stable orbits there except those tightly 
bound to one of the BHs. In principle, a small flux of particles to 
this region can produce the measured precession, but torques due 
to the anisotropy of the potential at larger distances are more likely 
to cause the precession. For comparison, general relativistic pre- 
cession of the orbit (not reproduced in the simulation) would be 
56° Myr -1 . We cannot yet say whether this precession compro- 
mises resonances that possibly induce the eccentricity growth, but 
will refer to this point in a future paper. 



6 DYNAMICAL FRICTION 

The orbital decay of a massive object within a galaxy down to its 
centre is well approximated by Chandrasekhar's dynamical fric- 
tion formula (Chandrasekhar 1943). However, the assumption of a 
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Figure 8. The semi-major axes a and eccentricities e at as functions of time, 
for three 10:1 binary simulations with initial eccentricity of 0.2. The dotted 
green lines are two different realizations with 2 x 10 5 particles while the 
solid blue line is a run with 5 x 10 6 particles. The first two simulations are 
arbitrarily terminated at t = 10000 while the last is stopped only when the 
eccentricity reaches 0.99, this time for a 10 8 M ., primary is 22 Myr from 
the beginning of the simulation. For this mass scaling, the timescale for 
coalescence due to gravitational waves at the end of the simulation is just 1 
Myr. 



uniform background in his classical treatment does not hold in real 
galaxies. Thus, corrections to the Coulomb logarithm are necessary 
to account for the changing background with radius. Just & Penar- 
rubia (2005) performed a detailed theoretical investigation of Chan- 
drasekhar's formula in the presence of a density gradient and gave 
an improved analytical formula for the Coulomb logarithm. Just et 
al. (201 1) took into account also self consistent velocity distribution 
functions, and made a comprehensive examination of the applica- 
bility of the new formula to sinking massive objects in a number 
of galaxy models, using high resolution A'-body and particle-mesh 
codes. Their results suggest a delay in the orbital decay with respect 
to the standard formula, which quantitatively varies according to 
the studied case. 

They notably give an explicit solution for the decay of a massive 
object moving on a 'circular' orbit. Their formula (equation 25 in 
their paper) is very general and holds for an arbitrary power law 
density profile (it does not hold for very flat cores where fast mov- 
ing stars contribute to most of the frictional drag; see Antonini & 
Merritt 2012). Here we bring their formula in our model's units and 
adjusted the parameters for an isothermal sphere; the radial evolu- 
tion is given implicitly by: 



: 52.918 x {Ei [21n (±/? )] -Ei [21n(|r)] } 



(27) 



The prefactor is an exactly calculable number. The special function 
Ei(jc) is called the exponential integral, it has real values only for 
x > 0. Thus, the domain of definition of equation (27) is r > 8. 
However, the assumption of a nearly circular orbit breaks well 
above r = 8. The angle between the velocity vector and the tan- 
gent can be derived by finding the radial velocity from equation 
(27): 



9 = 195.99° 



1 



■ln( 



(28) 



where this approximation hold only for small angles or large r. 

It is generally difficult to simulate a full infall of a compact object 
into a galaxy centre due to the collisional nature of the interactions 
and the low number density of particles which can generally be 
obtained at the outskirts of galaxy models. Just et al. (2011) used 
the particle-mesh code SUPERBOX (Fellhauer et al. 2000), which 
is collisionless and uses fixed time steps; This type of code, unlike 



Figure 9. This simulation of an equal mass binary initiated from Rq = 500 
(solid blue line) focuses on the dynamical friction regime. The theoretical 
curve (dashed black line) is the analytical formula of Just et al. (201 1) given 
by equation (27). The orbital decay is initially very well described by the 
formula, and the deviations at r < 300 are possibly due to inconsistency of 
the phase space distribution in the actual simulation versus the assumption 
of isothermal sphere used to derive the formula 



direct A'-body codes, allows a large particle number to be simulated 
in a relatively short time. Using the code we designed for BBHs, we 
can also study the early part of the inspiral, which is dominated by 
dynamical friction. Since our code resolves collisions between the 
BHs and stars, this study is complimentary to that of Just et al.. 

In Figure 9 we show a simulation of an equal mass binary ini- 
tiated from Rq — 500 (solid blue line) and the theoretical curve 
(dashed black line), equation (27). The initial velocity is the local 
circular velocity, but the initial velocity vector is tilted by 1.62°; 
this angle is obtained be substituting r = Rq in equation (28). If the 
initial velocity is purely tangential, then the spiral becomes slightly 
'eccentric'. A second run with a different realization (not shown) 
gave very similar results, including the position of the wiggles. It 
should be pointed out that this simulation is just a proof of concept; 
we do not expect the stellar distribution to be spherically symmetric 
at the early stages of a major merger. 

The orbital decay in our simulation is initially very well de- 
scribed by equation (27); deviations become significant below r ~ 
300, where mutual gravity of the two BHs is still negligible com- 
pared to the gravity of the bulge. These deviations are possibly due 
to the fact that the phase space distribution in the actual simula- 
tion (at the time and radius where the deviations occur) is no longer 
consistent with the assumption of isothermal sphere used to derive 
equation (27), in particular the velocity distribution might not be 
described well by a Gaussian. 



7 SUMMARY 

Using a conservation-based scheme, we were able to follow the 
evolution of a BBH from a wide separation (enclosed stellar mass 
greater than the combined BH mass) down to sub-parsec scale. Our 
code resolves star-BH collisions and can run with N > 10 6 stars on 
a desktop computer. We verified that our scheme yields convergent 
results which are independent of the number of particles, and the 
initial and boundary conditions. By performing scattering experi- 
ments on the inspiraling BBH, we will be able to extend Meiron 
& Laor (2010) and calculate the signature of the inspiral on the 
background stellar phase space distribution as a function of pro- 
jected position. This calculation improves on A'-body simulation 
by reducing statistical fluctuations and having no spurious relax- 
ation (and thus no loss cone refilling). We performed calculations 
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for both an equal mass binary and a 10:1 mass ratio. Our calcula- 
tions reveal: 

(i) The inspiral from a radius scale of tens of parsecs to the hard 
binary radius occurs on a time scale of a few million years for a 
10 8 Mq primary, with only a weak dependence of the timescale on 
the mass (°= Mg 29 , equation 19). 

(ii) The inspiral ends at a radius which is ~ 30 per cent smaller 
than the simple analytical estimate for the hard binary radius, and 
consistent with Merritt (2006). 

(iii) An equal mass binary inspiral leads to a nearly circular final 
orbit, regardless of the initial eccentricity. 

(iv) Eccentricity increases and coalescence due to gravitational 
wave emission will occur for a binary with a mass ratio of 10:1 in 
less than 10 8 years (xMg ,s ). If the stellar distribution is triaxial or 
rotating the lifetime of such systems is potentially shorter. 

While we used a static, spherically symmetric background poten- 
tial to account for star-star interactions, it is straightforward to ex- 
tended this method to treat more complicated cases such as an adi- 
abatically evolving potential (e.g. due to core depletion during the 
BBH inspiral), non symmetric stellar models and perturbers within 
the scattering method. 

This method can also be used to explore the time evolution of 
the statistical properties of the scattering events, and the exten- 
sion of the dynamical friction formulation to the hard binary stage. 
This study will help understand the mechanism which leads to the 
growth of eccentricity for an q ^ 1 BBH, and its decay for the q = 1 
case. 
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